A STABILIZED NITSCHE FICTITIOUS DOMAIN METHOD FOR 
THE STOKES PROBLEM 
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Abstract. We develop a Nitsche fictitious domain method for the Stokes problem starting 
from a stabilized Galerkin finite element method with low order elements for both the velocity 
and the pressure. By introducing additional penalty terms for the jumps in the normal velocity 
and pressure gradients in the vicinity of the boundary, we show that the method is inf-sup stable. 
As a consequence, optimal order a priori error estimates are established. Moreover, the condition 
number of the resulting stiffness matrix is shown to be bounded independently of the location of the 
boundary. We discuss a general, flexible and freely available implementation of the method in three 
spatial dimensions and present numerical examples supporting the theoretical results. 
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1. Introduction. A frequently encountered problem in practical applications 
of the finite element method is the generation of a high quality mesh conforming to 
the computational domain. For instance, the simulation of flow around an object 
embedded in a channel typically requires a mesh discretizing the domain surrounding 
the object. If the domain is complex, the mesh generation problem is highly non- 
trivial. Furthermore, the mesh must be modified or regenerated each time the object 
is translated, scaled or rotated, for example to study the lift or drag for different 
angles of attack. 

In fictitious domain finite element methods [HI [17l [22l [37] , the computational 
domain is instead represented by a, possibly regular, background mesh and an interior 
surface; this situation is illustrated in Figure |1.1| The mesh generation problem is 
thus essentially avoided. However, new challenges are introduced. The interior surface 
must be represented and the intersection of the surface and the underlying mesh 
computed, which is a complex task for three-dimensional domains. Moreover, the 
finite clement formulation, and hence also its analysis and implementation, involves 
elements of non-regular shapes induced by this intersection. 

In this work, we consider a Nitsche fictitious domain method for the Stokes prob- 
lem: find the velocity u : il C M. d — > R d and the pressure p : il — > R such that 



-Au + \7p = f mil, (1.1a) 
V-m = mil, (1.1b) 
u = g on r, (l-lc) 

where il denotes a bounded domain in R d , d = 2 or 3, with Lipschitz boundary 
r, and where / 6 L 2 (il) is a given body force and g € i? 1 / 2 (r) is a prescribed 



boundary velocity. To satisfy (1.1b), we assume that J T n ■ gds = where n denotes 
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Fig. 1.1. The stabilized Nitsche fictitious domain method presented in this work allows the 
simulation of Stokes flow around a possibly complex object (in this simplified illustration, a two- 
dimensional airfoil) embedded in a fixed background mesh. The object is defined by its boundary Y, 
and the computational mesh (here the cut mesh surrounding the airfoil) is defined as the intersection 
of the fixed background mesh and the outside (or inside) of the boundary Y . 

the outward pointing boundary normal. Moreover, we assume that J^pdx = to 
uniquely determine p. 

The fictitious domain method introduced in this paper is based on a least squares 
stabilized finite element method with low order finite element spaces. In particular, we 
consider both the case of continuous piecewise linear vector fields for the velocity and 
continuous piecewise linears for the pressure, and the case of continuous piecewise 
linear vector fields for the velocity and piecewise constants for the pressure. We 
prove stability and optimal a priori error estimates as well as optimal estimates for 
the condition number. These results rely on the introduction of stabilization terms 
for the jump in the normal gradients at faces associated with elements intersecting 
the boundary. Our method is closely related to a very recent report of Burman 
and Hansbo [T5J, but the analysis follows a different route. The present work also 
differs from that of Burman and Hansbo jl2j in that our methodology has been tested 
and implemented in three dimensions. Similar results have been obtained for elliptic 
boundary problems by Burman [10] , Burman and Hansbo and Johansson and 
Larson [55]. In a related work [3U], we present a stabilized Nitsche overlapping mesh 
method for the Stokes problem. 

A central and unique contribution of the current work is the full and general 
treatment of domains represented by arbitrary boundary triangulations embedded 
in three-dimensional tetrahedral meshes. This requires integration over arbitrary 
polyhedral domains resulting from the subtraction of the embedded domain from 
the background mesh. The intersection of the boundary and the background mesh 
is computed efficiently using techniques from computational geometry. The freely 
available implementation is based on, but extends that of, our previous work |29j . 

The remainder of this paper is organized as follows. In Section [2j we summarize 
the notation and assumptions used throughout this work. The novel Nitsche fictitious 
domain finite element formulation for the Stokes problem is then introduced in Sec- 
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Fig. 2.1. (Left) The computational domain Q is defined as the inside or outside of a given 
boundary F imposed on a fixed background meshl~*. (Right) The fictitious domain Q* is the union 
of the minimal subset T* C T* covering Q. 



tion [3j while Sections [4|{6] are devoted to its a priori error analysis. We prove that 
the condition number is bounded independently of the location of the boundary in 
Section [7J A brief summary of key implementation aspects is provided in Section [8j 
along with numerical investigations corroborating the theoretical results and an ex- 
ample demonstrating the applicability of the developed framework to complex 3D 
geometries. Finally, we provide some concluding remarks in Section [9] 

2. Preliminaries. The Nitsche fictitious domain finite element formulation in- 
volves integration over various geometric entities. We here define these entities and 
summarize the notation that will be used throughout this paper for computational 
domains, meshes, function spaces and norms. 

2.1. Computational domain and meshes. Let SI be an open, bounded do- 
main in M. d (d = 2,3) with Lipschitz boundary T. We assume that SI is a subset of 
a larger polygonal domain ST; that is, SI C fi*. We will refer to ST as the fictitious 
domain. Let T* be a shape-regular tessellation of SI* such that T D SI ^ for all 
T G T* ■ The mesh T* might be constructed fro m a la rger and easy-to-generate mesh 



T* by extracting a suitable submesh, cf. Figure 2.1 A facet F; that is, an edge in 
two dimensions or a face in three dimensions, of the mesh T* is labeled an exterior 
facet if it belongs to one element only (and is thus a part of the boundary of SI*) or 
an interior facet if it is shared by two elements. In the latter case, we denote the two 
elements shared by the facet F by TJ and Tp . The set of all exterior facets defines 
the boundary mesh d e T* , while the set of all interior facets defines the skeleton mesh 
diT*. 

Given T*, we may define the cut mesh T on SI as follows: 

T = {THO : T G T*}. (2.1) 

The corresponding boundary and skeleton meshes are defined accordingly by d e T — 
{F n H : F G d e T*} and d{T = {F n O : F G d,T*}. Note that T, d e T and diT 
consist of both standard (simplicial) elements and facets, and non-standard elements 
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/ \ r 

Fig. 2.2. The boundary zone of the fictitious domain. (Left) The background meshT* together 
with the cut mesh T. Observe that for the element associated with the node x, only a small fraction 
resides inside the domain Q. (Right) The elements in yellow are intersected by the boundary and 
therefore part of the mesh ■ Interior facets of elements intersected by the boundary (J 7 ^) are 
marked in green. 



and facets. We will occasionally refer to the former set as non-cut elements or facets, 
and the latter set as cut elements or facets. 

Next, let Tr* be the subset of elements in T* that intersect the boundary T: 

Tr* = {T e T* : T n T ^ 0} (2.2) 

and introduce the notation T^ for the set of all interior facets belonging to elements 
intersected by the boundary T: 

Fr={Fe diT* : T+ n T ^ V T F n T ^ 0}. (2.3) 

Figure [2~2| illustrates this notation. 

We assume that T* and the boundary T satisfy the following geometric conditions: 

• Gl: The intersection between T and a facet F € diT* is simply connected; 
that is, r does not cross an interior facet multiple times. 

• G2: For each element T intersected by T, there exists a plane St and a 
piecewise smooth parametrization <f> : St H T —> T n T. 

• G3: We assume that there is an integer N > such that for each element 
T G 7f* there exists an element V G T* \ Tr* and at most TV elements {T}f =1 
such that T x = T, T N = T and T n T i+1 e diT*, i = 1, . . .N - 1. In 
other words, the number of facets to be crossed in order to "walk" from a cut 
element T to a non-cut element T" C f2 is bounded. 

Similar assumptions were made by Burman and Hansbo [llj , Hansbo and Hansbo [18] 
for the two dimensional case and ensure that T is reasonably resolved by T* ■ 

2.2. Finite element spaces. We let the discrete velocity space Vu be the space 
of continuous, piecewise linear M d -valued vector fields defined relative to a specified 
mesh, and let the pressure space Qh consist of either piecewise constant or continuous 
piecewise linear elements, denoted by P ^' dc and P^, respectively. 

Here and below, let || • \\ s .n and | ■ \ s ,n denote the standard Sobolev norms and 
semi-norms on a domain Q for s G N. The corresponding inner products are denoted 
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by (v)s,n- For s — 0, the subscript s is omitted. Furthermore, we introduce the 
following mesh-dependent norms for the velocity: 

IIMII 2 = ||Vw|ft + ||^ 1/2 -y||r + \\h 1/2 n ■ Vv\\l, (2.4) 
IIHIl2 = l|V«||^ + ||/i- 1/2 «||^ (2.5) 



for the pressure: 



and for the product space: 



MI| 2 =IMI n + ll^ 1/2 <zllr, (2-6) 
Mil* = 119110.= (2-7) 



|||(^5)||| 2 = ||H|| 2 + |||g||| 2 , (2.8) 
\\\(v,q)\\\l = \\\v\\\l + \\\q\\\l (2.9) 

Note that the ||| • |||*-norms are defined on the fictitious domain Q* and therefore rep- 
resent proper norms for the discrete finite element functions. When mesh-dependent 
norms are applied to non-finite element functions on a domain f2, we always mean the 
evaluation of the norm on a tessellation T of 17. 

3. Finite element formulation. Before we present the Nitsche fictitious do- 
main method, we review a pair of well-established stabilized finite element formula- 
tions for the Stokes problem. These formulations are then extended to a Nitsche-based 
fictitious domain method. 

3.1. Stabilized Stokes elements. Let Vh and Qh be the velocity and pressure 
spaces introduced in the previous section defined relative to a standard conforming 
tessellation T of f2 and recall that Qh is defined to be either P^ or P°' dc . It is well- 
known that the mixed spaces Vh, X and Vh X P°' dc violate the inf-sup condition for 
the [£fQ(fi)] d x L 2 (tt)/M. variational formulation of the Stokes problem ( 1.1 1, and thus, 
are not stable in the Babuska-Brezzi sense [S]. Different strategies can be employed 
to compensate for the lack of stability [7J [TH1 HOI 123] > whereof consistently stabilized 
methods are among the most prominent [51 114). Here, we consider consistently stabi- 
lized discrete variational formulations of ( 1.1 1, with g = 0, of the following form: find 
(u h ,Ph) € V h x Q h such that 

A h (u hl p h ;v h ,qh) = Lh(v h ,qh) y(v h ,qh) € V h x Q h , (3.1) 
where the bilinear and linear forms Ah and Lh are defined by 

Ah(u h ,Ph;v h ,qh) = a h (u h ,v h ) + b h {u h ,q h ) + b h {v h ,Ph) ~ c h (u h ,Ph;qh), 

L h (v h ,q h ) = {f,v h )-$ h (q h ). (3.2) 

Here, ah and bh are the standard forms 

a h (u h ,v h ) = (Vuh,Vv h )n, (3.3) 

bh(vh,Ph) = -(y-Vh,Ph)n- (3-4) 
The stabilization form Ch is given by 

jO.dc 
h 

Pi Etgt hU-^u h + V Ph , Vq h ) T if Qh = Pt 



r ( . ^_ I A) Ef69,r h F {\ph], [qh\)F if Qh = P h , /„ 

c h (u h ,p h ,q h )-{ h 2,_ A „.. „„.n ;f n._pi ^ 
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where hx denotes the diameter of element T, hp denotes the average of the diameters 
of the elements sharing a facet F, [v] = v + — v~ is the jump in a function v over each 
facet F: v ± {x) = lim t _> + v(x±tn) for x £ F, and (3q and Pi are positive stabilization 
constants. In the case Qh = Phi this stabilization is also known as the pressure- 
Poisson stabilized Galerkin method. Note that -Am/j vanishes if Uh is piecewise 
linear and is only included to clarify that the method is indeed consistent. We will 
therefore simply write Ch(ph, qh) when only finite element functions are involved. The 
form <j>/j in (3.2) is, to ensure consistency, defined to be 

[Pi J2t€T h T\fi Vg h ) T if Q h =Ph- 

Since [qh] — for qh <G P^ and V((?^|t) = for qh € -P"' dc , we may express the two 
cases in a more compact notation: 

ch(ph,qh) = A) X! ^(bJ" + A E 4(%> v %)t, (3.7) 

Fed.T TeT 

*kW=ftE^(/.%^ (3-8) 
TeT 

3.2. A stabilized Nitsche fictitious domain method. Prior to stating the 
stabilized Nitsche fictitious domain formulation for the Stokes problem, we introduce 
the following forms with reference to the notation established in Section |2.1| 

a h (u h ,Vh) = (Vu h ,Vv h )n - (d n u h .v h )T - {d n v h , u h )r + "/(h~ 1 u h , v h )r, (3.9) 
bh(vh,Ph) = -(V ■ v h ,Ph)n + (n ■ v h ,Ph)r, (3.10) 
where d n v = n ■ Vu. Next, we introduce the velocity "ghost-penalty" form: 

i h {u h ,v h ) = /?2 ^2 h F ([d n u h ],[d n v h ]) F , (3.11) 

and the pressure "ghost-penalty" form: 

. , * fPoT, F &F* h F([ph],[<lh])F\a if Qh = P°' Ac , , qio v 

3h(J?h,qh) = < _ „ t.3na l r* n ■* ^. D i t 3 ' 12 ) 

(ALir 6 j ; h F ([d n p h ], [d n q h ])F if Qft = Pf\- 

Again, [v] = v + — v~ is the jump over each facet F, and n — rip is a fixed, but 
arbitrary, unit normal to the facet F. Here, A > and /?3 > denote additional 
penalty parameters. As before, we are allowed to rewrite (3.12) as a single form 
jh{Ph,qh) = jh,o(Ph, qh) +jh,i(Ph,qh) with j 0>h and j h .i denoting ( |3.12[ ) in the case of 
Qh = Ph' dc and Qh = P t l, respectively. 

We are now ready to state the Nitsche based fictitious domain method for the 
Stokes problem (LTJ. Let V h = V h {T*) and Q h = Qh(T*) be the finite element 
velocity and pressure spaces defined relative to T* . The variational problem reads: 
find (uh,Ph) e V h x Q h such that 

Ah(u h ,p h ;v h ,q h ) + Jh(uh,PhiVh,qh) =L h (v h ,q h ) y(v hl qh) € V/, x Q h , (3.13) 
where and Jh are defined by 

A h (u h ,Ph]V h ,qh) = a h (u h ,v h ) + b h {u h ,qh) + b h {v h ,Ph) ~ Ch(ph,qh), (3-14) 
Jh(u h ,Ph;v h ,qh) = ih(u hl v h ) - jh(ph,qh), (3.15) 
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where the forms Ch and $/j are defined as in (3.7) and (3.8) (relative to the cut mesh 
T). The form Lh is given by 



L h (v h ,q h ) = (f,v h ) n + (9,7/1 1 v h - d n v h +q h n) r - $ h (q h ). 



(3.16) 



Remark 3.1. The "ghost-penalty" defined in (3.11) was introduced by Burman 



and Hansbo Jltf to formulate a first- order convergent fictitious domain method for the 
Poisson problem. The role of the ghost-penalty is to augment the bilinear form ah by 
accounting for small elements |T D fi| <C |T|, T G T* in the vicinity of the boundary 

r. 

Remark 3.2. In the Stokes problem, the stabilization form c/j acting on the 
pressure also has to be augmented. Depending on the pressure discretization, this can 
be achieved in different ways. In the case of Qh = Ph' d °, a similar ghost-penalty was 
presented by Becker et al. [6] to propose a finite element method for incompressible 
elasticity problems with discontinuous modulus of elasticity. To motivate the ghost- 



penalty (3.12) when Qh — P^, one may consider the stabilization terms fafi(Vp, V<z)t 
and h^ r (f 1 Vq)T as a locally scaled version of a Poisson equation and apply (3.11). 
In Lemma \5.1\ we will reveal the basic structure behind the augmentation terms and 
also present a generalization to higher-order elements. 

4. Approximation properties. Before we proceed with the a priori error anal- 



ysis of the method proposed in Section 3.2 we summarize here some notation and 
useful inequalities that will be used throughout Sections [5] and [6j In what follows, V^ 
and Vh denote some finite element spaces consisting of piecewise polynomial functions 
defined on T* and T respectively, but it should be clear that we have mainly Vh = Vh 
or Vh = Qh in mind. The constants C involved in the inequalities will only depend 
on ft or fi* , the regularity of the relevant function spaces, the shape-regularity of T* , 
and possibly the polynomial order of Vh] in particular, the constants C do not depend 
on h. 

4.1. Trace inequalities and inverse estimates. We recall the following trace 
inequalities for v € i/ 1 (£!*): 



IMIor < C(h- 1/2 \\v\\ T + 4 /a ||V«|| T ) VT G T, 
Httst < C{h T 1/2 \\v\\ T + 4 /2 I|Vw||t) VT e T. 



(4.1) 
(4.2) 



See Hansbo and Hansbo [TH] for a proof of (4.2). We will also need the following 
well-known inverse estimates for Vh G Vh' 



\\Vv h \\ T < CTi^lkHr VTeT*, 
Ih^n-VvhWF^CWVvhWr VTgT*, 



(4.3) 
(4.4) 



For proofs, we refer to Quarteroni 
the boundary parts T (~\T: 



Moreover, we will need a version of (4.4) for 



|/i 1/2 n-V^||rnT<C||V^||T VTgT*, 



(4.5) 



which was proved under assumptions similar to Gl - G2 by Hansbo and Hansbo [18 . 
We note that for f/i € 14, ft 6 Qh, we have the two estimates: 



kill «S G\\\v h \ 



Mil *s c\\\qh\ 



(4.6) 



which can easily be deduced by (4.5), and by combining ( |4.2| and (4.3) 
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4.2. Interpolation estimates. In order to construct an interpolation operator 
L 2 (fl) Vh, we recall that there is a linear extension operator £ : H s (ft) —> H s (fl*), 
s ^ 0, such that 

||£u||.,n. < C||u||., n . (4.7) 

See Stein [33J for further details. Let -K* h : L 2 (Vt*) — > V*_ be the standard Scott-Zhang 
interpolation operator |32| and recall the interpolation error estimates 

\\v-Tr* h v\\ r , T ^Ch s - r \v\ sMT) , 0<r<s<2 VT £ T, (4.8) 
\\v-Tr* h v\\ rtF ^Ch s - r -^ 2 \v\ sMT} , 0^r^s^2 VFeftT*, (4.9) 

where uj(T) is the patch of neighbors of element T; that is, the domain consisting of 
all elements sharing a vertex with T. Next, we define 717, : L 2 (fl) —> Vt as follows: 

n h v = ir* h £v. (4.10) 



Note that n^v is now defined on 17*, and in particular on 12 C 17* 



The stability estimate (4.7) together with the interpolation error estimates (4.8) 



and (4.9) for the Scott-Zhang interpolation operator imply the following interpolation 



estimates: 



\\v-n h v\\r,T^Ch s - r \v\ sMT) , 0<r<s<2 VTeT, (4.11) 

\\v-n h v\\ r ,F < Ch s - r -^ 2 \v\ StUl{T) , 0<r< S <2 VFeftT. (4.12) 



We now return to our specific finite elements spaces Vh and Qh- For the energy 
norm, we have the following interpolation error estimates: 



Lemma 4.1. For the interpolation operator n h defined by (4.10), there is a con- 
stant C > such that for all v £ [H 2 (fl)] d and all q £ i? 1 (SI): 



\\\v - ir h v\\\ < Ch\v\ 2 ,n, 
(v - n h v, q - TT h q)\\\ < Ch(\v\ 2 ,n + 



(4.13) 
(4.14) 



Proof. We only sketch the proof. First use the trace inequality (4.2) to esti- 
mate the boundary contributions in terms of element contributions. Then apply the 
interpolation error estimate (4.8), and finally the stability estimate (4.7). □ 

In addition to the interpolation estimates, we will need the following continuity 
property of the extended interpolation operator with respect to different norms: 

Lemma 4.2. Assume v £ [H^(Q ,)] d and let ix h : \H 1 {Vl)\ d -s- V h (T*) be the 
interpolation operator defined in (4.10). Then there is a constant C > such that 



Whv\\\* < C\\v\\ 



i,n- 



(4.15) 



Proof. By definition we have |||7ivu|||* = ||V7ivu 



h-^TVhvWf. The bound 



for the first term on the right-hand side follows immediately by the boundedness of 
and the continuity of the extension operator £. To estimate the second term, we use 
the fact that £v\y = for v £ [HQ(fl)] d , the trace inequality (4.2), the interpolation 
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h 



Fig. 5.1. Controlling the L 2 -norm \\v\\t„ of a finite element function v on a barely intersected, 
'fictitious" element To by ||^||t3 md boundary zone jump-penalties. Starting from To, each term 
can be estimated by the neighboring term \\v\\ T 2 when a sum of jump-terms of the form 

2J+ 1^|I? 



i+i 1 



is added. 



estimate (4.8 1 and continuity of £ again: 



IrnT 



\h 1,2 TT h v\\l= h T l hhv\\l nT = Y h T 1 \\n h v-£v\ 
TeT* TeT* 

< ^ ^t 1 {^ l T 1 \\ 7r h v — £v\\t + ^t||V (ir h v — £v) Hy) 

TeT* 

<C\\£v\\l n ,^C\\v\\l n . 



5. Stability estimates. In this section, we demonstrate that the bilinear form 



defining the stabilized Nitsche fictitious domain variational formulation (3.13) indeed 
satisfies the inf-sup stability condition in the Babuska-Brezzi sense. 



5.1. The role of the boundary zone jump-penalties. As a first step, we 



show how the jump-penalties (3.11) and (3.12) in the boundary zone J-"p contribute 
to control the norms of and p^ on the entire fictitious domain T* . We start with 
the following lemma. 

Lemma 5.1. LetT = {T} be a tessellation consisting of shape-regular elements T 
and let T\, Ti € 1" be two elements sharing a common face F. Furthermore, let v be a 
piecewise polynomial function defined relative to the macro-element T = 7\ U T 2 . Let 
Vi be the restriction of v to Ti for i = 1,2. Then there is a constant C > 0, depending 
only on the shape-regularity ofT and the polynomial order p = max(ord(«i), ord^)) 
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of v, such that 



\H 2 Tl <C \ \\v\\ 2 T2 +J2h 2j+ \[dLv],[div])F) , (5.1) 

where d^v = J2\ a \=j D a v(x)n a for multi-index a = (ax, . . . , ad), \a\ — J2i a i an d 
n a =n^n% 2 ■■■n a d d . 

Proof. For a given point x G T±, we write xp = xp(x) for the normal projection 
of x onto the plane defined by the face F. Note that the area \T\ p\ of all projected 
points in T\ is bounded by |F| up to a constant by the shape-regularity assumption. 
For i — 1,2, and since ord(ui) p, we may express (the extensions to T of) Vi in 
terms of its Taylor-expansion around xp: 

Vi{x)= y . {\x-xp\n) a , 

where n is the unit normal vector of F pointing towards T±. Subtracting the two 
Taylor expansions, we find that 

Vi(x) = v 2 (x) + } j (\x-x F \n) a . 

Next, integrating over T\ with respect to x, taking squares and applying the Cauchy- 
Schwarz inequality yield 

'"j II v, ( ' ( E / ([D a v(xp(x))]n a ) 2 h 2 ^ dx\ , 

with h the maximal element diameter. From the assumption of shape regularity, a 
change of variables, and the definition of d° n , it follows that 



Mil, < C | \\v 2 f Tx +E / [«y)] 2 ^ +1 dy 



Finally, as the two norms ||w2||ti an( i ||^2||t 2 are equivalent, again by shape regularity, 



we obtain the desired inequality (5.1). □ 



Remark 5.1. The previous lemma is a key observation for proving stability and 



a priori error estimates for the fictitious domain formulation (3.13) as it lays the 
foundation for how to control certain norms on the fictitious domain £1* in terms of 
norms computed only on f2 and appropriate jump-penalties in the intersection zone 
J-~p . We are now in a position to state the following proposition: 

PROPOSITION 5.1. Let £1, fl* and F£ be defined as in Section [Q} There is a 
constant C > such that the following estimates hold. 
Forallv h eV h (T*): 

\\Vv h \\l, sc C(\\Vv h \\l + £ hp([d n v h ], [d n v h ]) F ) < C\\Vv h \\ 2 n , (5.2) 
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and for all qt € Ph' dc (T~ 



IMIo» < C(\\q h \\l + h F ([n-q h ],[n-q h ]) F ) < C|| 



l|2 



while, for all g% g -P/j(T*) 



Ikhlln- < C(ll%lln + E ^f([M,M) f ) 



(5.3) 



(5.4) 



Proof. We start with the first inequality of (5.2) 
into sums over non-cut and cut elements 



Decompose the norm over f2* 



Let To € Tp* be a cut element. By the 
geometric condition G3 (cf. Section 2.1), there exists a Tjy C f2 and at most N — 1 
elements T; £ Tp* and facets T,_i fl Tj = G J-p that have to be crossed in order 
to traverse from To to T/v- By the shape-regularity of the mesh, each facet F g .Fp 
will only be involved in a finite number of such crossings. Applying Lemma |5.1| with 
each component of Vi^ as v, iteratively to each neighboring pair {Ti_i, Tj} yields the 
desired estimate. 

Th e firs t inequalities of (5.3) and (5.4) follow by the analogous argument: apply 
Lemma 



5.1 



to qh and recall that [qh\F — for qh e P^. 



The second inequalities of ( 5.2 H 5.4) rely on the shape regularity, allowing us to 



bound hp by h, and the trace and inverse estimates of Section |4~T| applied to each 
facet of the boundary zone sums. The upper bounds immediately follow. □ 

Remark 5.2. Burman and Hansbo Ulf presented the analogous result to (5.2 ) for 
the Poisson problem with continuous piecewise linear finite elements. The formulation 
given here, together with Lemma \5.1\ reveals the basic structure of jump-penalty-based 
stabilization terms for fictitious domain formulations and can be applied to various 
types of norms and elements, including higher-order elements. 

5.2. Stability estimates and the inf-sup condition. The main result of this 
section, Theorem |5.1[ is the inf-sup stability of the bilinear form Ah + Jh , occurring 



in the stabilized Nitsche fictitious domain variational formulation (3.13 1, with respect 
to the III 



norm (2.9). 



We begin by establishing the properties of the separate contributions to the bi- 



linear form. First, the form cf. (3.9) augmented by ih cf. (3.11) is continuous and 
coercive with respect to the norms ||| • ||| and || 
constants c > and C > such that 



|llj . More precisely, there are 



a h (v,w) sc C\\\v\\\ \\\w\\\ \/v,we [H 1 {T)] d , 
a h {v h ,w h ) + i h (v h ,w h ) ^ C\\\v h \\\* \\\wh\\\* Vv h ,w h G V h (T*), 

c\\\w h \\\l sC a h (w h ,w h ) + i h (w h ,w h ) Vw h £ V h {T*). 

Next, we show that bh is continuous with respect to the relevant norms. 



(5.5) 
(5.6) 
(5.7) 



Proposition 5.2. Let b h be defined by (3.10). There is a constant C > such 



that 



b h (v,q)^C\\\v\\\ 
b h {v hl q h ) ^ C\\\v h \ 



\1h\ 



V(M)e|ff Wxi 2 P), 
V(v h ,q h ) S V h (T*) x Q h (T*)- 



(5.8) 
(5.9) 
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Proof. The bound (5.9) follows from the definitions of bh and the 



and a subsequent use of (4.2 1 and (4.3). □ 



norm, 



The next lemma gives a fictitious domain adapted version of a "bad inequal- 
ity" often used in Verfiirth's trick and in proofs for some classical, stabilized 
schemes [14]. 

Lemma 5.2. There are positive constants C\,Ci such that for each qh £ Qh{T) 
there exists a Vh € V^(T*) satisfying 



bh(v h ,q h ) 



\ V h\ 



>C 1 \\q h \\ a -C a ((Y i h^\\Vq h \\^ 2 + ( E MM!") 1 



/ 2 1 



(5.10) 



TeT 



Proof. Let qh € Qh(T) be given. There exists a v £ [H 1 (fl)] d and a constant 
Ci > such that divv = g h and C*i||f ||i,o < || divt>||o [15]. Map u i-> 7T/,w € V h (T*) 
by the extended interpolation operator cf. (4.10), and denote e h = Tr h v — v. It follows, 
using the definition of bh , that 



bn{^hV,qh) = b h (e h ,q h ) +b h (v,q h ) ^ b h (e h ,q h ) + Ci\\v\\ lt n\\qh\\a- 
Moreover, integrating by parts on each element TeT yields 
bh(eh,<lh) = -(dive h ,q h ) u + (n ■ e h ,q h ) r 

= E ( e '" V * i ) T + E {n- e h ,[q h }) F , 



(5.11) 



(5.12) 



TeT FediT 

while the Cauchy-Schwarz inequalities give 

Ek.%) t > -(E ^t 2 ii^iiI) 1/2 (E 4iiv%iii) i/2 , 

TeT TeT TeT 

E [©.])*■£-( E ^ii^Ht) 172 ! E Mik]iiD 1/2 - 

FediT Fed.T Fed % T 

Since h^We^r ^ C\\ y\\ lM (T) by (|4.11[ ) and fr F 1/ 2 ||e/ t || f < C||w|liMT) by ( |4.12[ ), we 
obtain by combining (5.12) with (5.13) and (5.14): 



(5.13) 
(5.14) 



bh(e h ,q h ) 



>-&IMkn((£4l|vrf) 1/2 + ( E M 

V TeT Te9iT 



MIIt) 



(5.15) 



Finally, combining (5.11) with (5.15), and recalling that |||7r/i»|||» ^ C'||w||i,n by 
Lemma 4.2 yields (5.10) with Vh — ithV. □ 



Using the stability estimates for and bh, we may now prove the following inf-sup 
stability estimate for for Ah + Jh- 

Theorem 5.1. There is a constant ca > such that for all (v,h,Ph) <= Vh x Qh = 
Vh{T*) x Qh(T*): 



A h (uh,Ph;v h ,q h ) + Jh(uh,Ph,;Vh,qh) ^ s 

SUP mT/ ViTi ^ CA|||(«fc,P/,) 

(»fc,9h)e%xQ fc \\\(Vh,qh)\\\* 



(5.16) 



Proof. The proof of Theorem 5.1 follows the proof by Franca et al. [14] . using the 
appropriate norms and Proposition 5.1 in combination with Lemma 5.2 Let (uh,Ph) 
be given. 
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First, choose (vh,qh) to be (— Wh,0) where —Wh satisfies (5.10) for the given ph- 
In addition, scale Wh such that = ||p^||n- For the sake of readability, we write 



k T (Ph) = (J2 h T\\VPh\\ 2 T) 1/2 + ( E h F\\\Ph}\\%) 1/2 - 
TeT Fed,T 



(5.17) 



With this choice of test functions, applying (5.6| and (5.101, and Cauchy's inequality 
with e give 

(A h + J h ){u h ,p h ; -w h ,0) = -a h (u h ,w h ) - i h (u h ,w h ) + b h (w h ,-p h ) 
^ -ClllwfclH* \\ph\\n + Ci\\ph\\n ~ C 2 kr(ph) \\Ph\\n 



> -^- e \\\uh\\\l + (Ci-e(C + C 2 ))\\p h 



Note that by definition kf(ph) ^ Kch(j>h,Ph) for some positive constant K depend- 
ing on Po and j3\. In combination with choosing e such that (C\ — e(C + C2)) > 0, 
this gives 

(A h + J h )(u h , Ph ; -w h , 0) > -C*|||M h |||2 + CilKHn - C 2 c h (p h ,p h ). 



Second, we take test functions (vh,qh) = ( U h,—Ph) which, using (5.7), gives 

(A h + Jh)(uh,Ph;u h ,-ph) = a h (u h ,u h ) + i h (u hl u h ) + c h (ph,Ph) +3h(Ph,Ph) 

> c\\\u h \\\l + c h (ph,Ph) +jh(Ph,Ph)- 

In total, for any 8 > 0, we have 

(A h + J h )(u h ,p h ;u h ,-p h ) +S(A h + Jh)(uh,Ph\ ~Wh,0) 

> (c - (JCtyllufclll^ + (1 - SC 2 )c h (p h , Ph ) + 5C 1 \\p h f a + j h (jp h ,p h ). 



Moreover, by (5.3) and (5.4), there exists a positive constant D such that 
\\Ph\\n+3h(ph,Ph) > D\\ Ph \\ 2 n , =D\\\p h \\\l 



(5.18) 



Finally, we conclude that by a suitable choice of S > 0, there is a positive constant ca 
such that (vh,qh) = (uh — Sw^, —ph) satisfies 

{A h + J h )(u h ,p h \v,q h ) ^ Ca(|||«/i|||^ + IIKIII*), 

which proves the desired estimate. □ 

6. A priori error estimate. Before we formulate the main a priori error esti- 
mate, we state two lemmas about how the stabilization form Jh affects the Galerkin 
orthogonality and the consistency of the total form Ah + Jh- Let Vh — Vh(T*) and 
Qh = QhiT*) throughout this section. 

Lemma 6.1. (Weak Galerkin orthogonality). Let (u,p) e [H 2 (£l)] d x if 1 (12) be 
the solution of the Stokes problem (1.1) and let (uh,Ph) be the discrete solution of the 
corresponding stabilized Nitsche fictitious domain formulation (3.13). Then, 

A h (u-u h ,p-p h ;v h) q h )- J h (u h) p h ;v h) q h ) = V(v h) q h ) e V h X Q h . (6.1) 
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Proof. The identify follows immediately from the fact that the solution (u,p) 
satisfies A h (u,p; v h ,q h ) = L h {v h , q h ), as defined by ( |3.14| ) and ( |3.16| ), for all (v h , q h ) £ 
x Q h . □ 

The ghost penalty part in involves the evaluation of d n qh on facets and 



therefore the variational formulation (3.13) is per se not consistent with (1.1) since 



we only assume that q £ The next lemma shows that this consistency error 

will not affect the convergence order. 

Lemma 6.2. (Weak consistency) Assume that u £ [H 2 (£l)] d and p £ 1 (SI) and 
let TTh be the interpolation operator defined by (4.10). Then for all {vh,qh) € Vh x Qh 
it holds that 



J h (ir h u,ir h p;vh,qh) ^ Ch(\u\ 2 ,n + \p\i,n)\\\(v h ,q h ) 



(6.2) 



Proof. By definition (3.15) 



J h (TThU,Tr h p;v h ,qh) = i h {iThU,Vh) - jh,o(^hP,qh) - jh,i(^hP,qh)- 



By the continuity assumption on u, ih(£u, Vh) = 0. So, by the definition of iih (4.10), 
the inverse inequality (4.4) and the interpolation estimate (4.8), and the continuity 
of £, we obtain 

ih(ir h u,v h ) = ih{n* h £u - £u,v h ) < C7i|£u| 2 ,n* ||V«h||n» s? C/i|«|2,n|||«/i|||*- 
Similarly, by the continuity assumption on p; the trace inequality (|4.1[), the inverse 



estimate (4.3) and the interpolation estimate (4.9); and the continuity of £: 
jh,o( n hP>Qh) = jhfi( n h £ P- £ P>1h) < Ci|£p|i,r2*IMIr2* < C%|i,n||k?s| 



Finally, to estimate jhi, we use (4.4) and (|4.3|; the boundedness of the Scott 



Zhang interpolant, and the continuity of £ to obtain 

jh,i{^hP,qh) «S C , /i||V7r /l p|| *|klk2* < C'/i|p|i i n|||g/ l |||*. 

Combining the three estimates yields the result ( |6.2[ ). □ 

Theorem 6.1. (A priori error estimate) Let (u,p) £ [H 2 {Vt)] d x if 1 (ft) be the 



solution of the Stokes problem (1.1) and let (uh,Ph) be the discrete solution of the 



corresponding stabilized Nitsche fictitious domain formulation (3.13). Then, there is 
a constant C > such that 



[u 



u h ,p~ph)\\\ s? Ch(\u\ 2 ,n + \p\i,a) ■ 



(6.3) 



Proof. Clearly, by the triangle inequality and (4.6): 

|||(U - U h ,p-ph)\\\ ^ \\\{U - TT h U,p- TThp)\\\ + C\\\{ir h u - U h ,TThP-Ph) 



Lemma |4 . 1 1 provides the desired bound for the first term on the right-hand side above. 
It is therefore enough to show that the discrete error (it^u — Uh, ithP — Vh) satisfies 
the error bound in (6.3). 

By Theorem 5.1 there exists a {vh,Ph) such that = 1 and 



CA\\\(K h U - u h ,TT h p- Ph)\\\* 

< A h (u h - n h u,p h - n h p;Vh,qh) + Jh(u h - n h u,p h - ir h p;Vh,qh) 

= A h (u - 7r h u;p - TT h p;v h ,q h ) - J h {^hU^hP\ v h ,qh), 
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where the last equality follows by the weak Galerkin orthogonality (6.1 1. Recalling 
the definition of Ah, we may write 

A h (u - ir h u,p- ir h p; v h , q h ) = a h (u - ir h u, v h ) 

+ b h (u - TT h u,q h ) + b h (v h ,p- ir h p) 
+ c h {u- ir h u;p- n h p, q h ). 



We use the stability estimate (5.5) for ah and (5.. 
estimate the first three terms: 



for b h ; and (4.6) and (4.14) to 



a h (u - n h u,v h ) + b h (u - 7rvu, + b h (v h ,p - ir h p) 

sC C(\\\u - ir h u\\\ \\\vh\\\* + \\\u - TT h u\\\ \\\q h \\\* + \\\vh\\\* \\\p ~ ^hP\\\) 

<C , /»(|u| 2l n + lp|i,n)|||(wfc,g/.)|IU. 

Using the fact that Air^u = locally and applying the Cauchy-Schwarz inequality, 
we may estimate the remaining term Ch = Ch(u — nhu;p — irhP, qh) by 



Ch=Po ^2 hF(.\p-ThP],[qh])F 

FediT 



/?! 2^ h^(-Au + V(p - TThPh),Vqh)T 



TeT 



hi\\Au 



2 N 1/2 
\T 



TeT 



(E^HV(p-7r h p)|||) 1/2 )( 



C( h F \\[p - n hP ]\ 

Fed.T 



TeT 

1/2 



E 

TeT* 



1/2 



)( 



FediT* 



h F\\[qh}\ 



1/2 



^Ch(\u\ 2 ,n + \p\i,n)\\\(v h ,qh) 



where we used the trace inequality (|4.1|) and inverse estimate (|4.3|) for the last term 
to pass to H^Hr* 



\1h\ 



Collecting all terms and applying the weak consistency 

(6.4) 



estimate (6.2) for Jh, we conclude that 

111717,14 - u h \\\ + \\\iThP - Ph\\\ < Ch(\u\a,n + \p\i,n), 
since \\\(vh, qh)\\\* =1-0 



7. Condition number estimate. Following the approach of Ern and Guer- 
mond [13] . we now provide an estimate for the condition number of the stiffness 
matrix associated with the finite element formulation presented in Section |3.2| In 
particular, the estimate shows that the condition number is uniformly bounded by 
Ch~ 2 independently of the location of the boundary T relative to the background 
mesh T*. 

First, we introduce some basic notation including the definition of the condition 
number. Let {<fi}f = i be a basis for some finite element space Vh- Then the expansion 
J2iLi Vi^i defines an isomorphism C : Vh — >• R w such that Cvh — V, where 



V = [Vi . . . V N ] T . We let (V, W) N = Y^Li Wi denote the inner product in R N and 
\V\% = (V, V)n the corresponding Euclidean norm. 



We introduce the stiffness matrix A such that 

(AV,W) N = A h (v h ,qh;wh,r h ) + Jh{v h ,qh;Wh,r h ) 



(7.1) 



for all v h , w h € V h and all q h , r h £ Qh where V = C(v h , qh) and W = C(w h , r h ). Since 
we consider the Stokes problem for an enclosed flow with the velocity prescribed on 
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the entire boundary T, the solution is only determined up to a constant pressure 
mode. Consequently, the matrix A is singular with kernel ker(_4) = span{C(0, 1)}. 
Throughout the remaining part of this section, we therefore interpret A as the bijective 
linear mapping between the — > M. N , where R N denotes the quotient space M. N = 
~BL N / ker(.4) and R N = im(A) = ker(A) 1 - the image space (note that A is symmetric). 
The condition number is defined by 

k(A) = \A\ n \A- 1 \ n , (7.2) 

with the operator norm 

I /II \ AV \ N f>7* 

\A\ N = sup -r— — . (7.3) 
veI N \o \ v \ N 



Equivalcntly, the operator norm \A\ may be defined by 

\A\ N = sup sup f ■ (7.4) 

For a conforming, quasi-uniform mesh T with mesh size h and a finite element 
space Vh defined on T, it is well known that there are constants c^ > and C M > 
only depending on the uniformity parameters and the polynomial order of Vh such 
that the following equivalence holds: 

c»h d,2 \V\ N ^ \\v h \\ ^ C^h d / 2 \V\ N Vv h eV h . (7.5) 

The following two lemmas are concerned with an inverse estimate and a Poincare 
inequality for the appropriate norms. 

Lemma 7.1. There is a constant Ci > such that 

IIKIH* s? Ci/T^KIIq. V« h € V h , (7.6) 

WKv^qhM^Cjh^Uv^qh)^ V(v h ,q h ) e V h x Q h . (7.7) 



Proof. By definition = || Vt> h\\n* +\\h~ 1 < 2 v h \\r- Hence, the inequality ( 7.6 1 

follows from the applying the inverse estimate (4.3) to the first term and the trace 
inequality (4.2) and subsequently ( |4.3[ ) to the second term. The second estimate (7.7) 
is a simple consequence noting that 1 ^ Ch^ 1 diam(fi). □ 

Lemma 7.2. (Poincare inequality) There is a constant Cp > such that 



\vh\\n* «S CplllufclH* Vv h e V h . 



(7.8) 



Proof. First we observe that \\v h \\n* = \\vh\\n + \\vh\\n*\n «S \\vh\\n + \\vh\\r*- 
To estimate ||f^||r2, we apply a variant of the standard Poincare inequality, valid for 

\\vf Q ^C(\\yvf n + \\vf r ). 

Since ||f||r ^ C||' l ~ 1 ^ 2 ' u lir\ we conclude that ||uft||n ^ Clll^hlll*- Using this estimate 
and the definition of Hl^hHI*, a bound for the remaining term ||i>/i||T r * can be obtained 
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as in the proof for (5.4) (noting that hp ^ Ch F )- 

\\v h f T *^C{\\v h \\l+ hU[d n v h ],[d n v h }) F ) 

< C(\\h- l ' 2 v\\l + \\Vv h f n + h F ([d n v h ],[d n v h ]) F ) 



Fen 



The last two terms are bounded by HVi^Hn* by (5.2), thus yielding (7.8). □ 

Finally, we state the continuity of the overall form Ah + J/, with respect to the 
norm ||| • HI*: 

Lemma 7.3. There exists a constant C A such that for all (vh,qh), (^h,Th) G 
V h x Q h 



A h (v h ,q h ;w h ,r h ) + J h (v hl q h ;w h ,r h ) sC C A \\\(v h ,q h ) 



(7.9) 



Proof. Because of the continuity estimates (5.6) and (5.9), it only remains to 



estimate the contribution Ch(qh,fh), which follows the same lines as in the proof of 
Theorem 16. II □ 

We are now in the position to state the main result of this section. 

Theorem 7.1. The condition number of the stiffness matrix A associated with 



the Nitsche fictitious domain method ( 3. 131) satisfies the estimate 



k(A) ^ Ch~ 



(7.10) 



Proof. Recalling the definition of the condition number in ( 7.2 ), the proof consists 
of deriving estimates for |-4|at and | ^4. 1 1 . By definition, for all V, W £ M. N \ {0}, 



(AV, W) N = A h (v h ,q h ;w h ,r h ) + J h (v h ,q h ;w h ,r h ) 
^C A \\\(v h ,q h )\\U-\\\(w h ,r h )\\U 



^C A Cjh- 2 \\{v h ,q h )\\ n ,-\ 

< c A c 2 c 2 h d - 2 \v\ N \w\ N , 



where the inequalities follow from the continuity of A/ l + Jh, the inverse estimate (7.7) 



and finally (7.5). Thus 



|^|jv < c A c]clh d - 2 . 

Similarly, for all V <G R N \ {0}, there exists a W such that 

(AV, W) N = A h (v h ,q h ;w h ,r h ) + Jh(v h ,qh;w h ,r h ) 
>c A \\\(v h ,q h )\\U-\\\(w h ,r h )\\\* 
^ Cp 2 c A \\(v h ,q h )\\ n . ■ \\(w h ,r h )\\n* 
^ C p 2 c A c 2 h d \V\ N \W\ N , 



(7.11) 



where the inequalities follow from the inf-sup estimate (5.16), the Poincare inequal- 



ity (7.8), and finally Moreover, 



\av\ n = su. p (av,y)\y\^ ^ (Avwwij/ 1 > c\v\ N , 

Y 
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where C = C P 2 CAC 2 l h d . Letting V = A 1 Y, which is allowed since A 1 is indeed 
invertible on the reduced space, and rearranging the inequality, we obtain |.A _1 3^|./v ^ 
C^IYU for all Y, and so 



IA^In s? C 2 P c^c- 2 h- 



(7.12) 



Combining (7.111 and (7.12), we obtain the desired estimate 



C C 2 

k(A) = \A\ N IA^In < CjC 2 P — ^h- 2 



CA Cf, 



8. Numerical examples. 

8.1. Software for fictitious domain variational formulations. The assem- 
bly of finite element tensors corresponding to standard variational formulations on 



conforming, simplicial meshes, such as (3.1), involves integration over elements and 



possibly, interior and exterior facets. In contrast, the assembly of variational forms 



defined over fictitious domains, such as (3.14), (3.15) and (3.16), additionally requires 



integration over cut elements and cut facets. These mesh entities are of polyhedral, 
but otherwise arbitrary, shape. As a result, the assembly process is highly non-trivial 
in practice and requires additional geometry related preprocessing, which is challeng- 
ing in particular for three-dimensional meshes. 

As part of this work, the technology required for the automated assembly of 
general variational forms defined over fictitious domains has been implemented as part 
of the software library DOLFIN-OLM. This library builds on the core components 
of the FEniCS Project [351 [57], in particular DOLFIN [35], and the computational 
geometry libraries CGAL [T] and GTS [2J. DOLFIN-OLM is open source and freely 



available from http://launchpad.net/dolfin-olm 



There are two main challenges involved in the implementation: the computa- 
tional geometry and the integration of finite element variational forms on cut cells 
and facets. The former involves establishing a sufficient topological and geometric 
description of the fictitious domain for the subsequent assembly process. To this end, 
DOLFIN-OLM provides functionality for finding and computing the intersections of 
triangulated surfaces with arbitrary simplicial background meshes in three spatial 
dimensions; this functionality relies on the computational geometry libraries CGAL 
and GTS. These features generate topological and geometric descriptions of the cut 
elements and facets. Based on this information, quadrature rules for the integration 
of fields defined over these geometrical entities are produced. The computational ge- 
ometry aspect of this work extends, but shares many of the features of, the previous 
work [21], and is described in more detail in the aforementioned reference. 

Further, by extending some of the core components of the FEniCS Project, in 
particular FFC [24] [28] and UFC [4] , this work also provides a finite element form 
compiler for variational forms defined over fictitious domains. Given a high-level 
description of the variational formulation, low-level C++ code can be automatically 
generated for the evaluation of the cut element, cut facet and surface integrals, in 
addition to the evaluation of integrals over the standard (non-cut) mesh entities. The 
generated code takes as input appropriate quadrature points and weights for each cut 
element or facet; these are precisely those provided by the DOLFIN-OLM library. 

As a result, one may specify variational forms defined over finite element spaces 
on fictitious domains in high-level UFL notation [3], define the background mesh T* 
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and give a description of the surface T, and then invoke the functionality provided 
by the DOLFIN-OLM library to automatically assemble the corresponding stiffness 
matrix. In particular, the numerical experiments presented below, corresponding to 



the variational formulation defined by (3.14), (3.151 and (3.16), have been carried out 
using this technology. 



8.2. Convergence rates. To corroborate the theoretical error estimate (6.3) by 
numerical results, we consider a basic test case with a manufactured exact solution 
and compute the errors in the velocity and the pressure approximations on sequences 
of refined meshes. To this end, let fl = [0, l] 3 with Y — <9f2. To examine the conver- 
gence of the Nitsche fictitious domain method, we apply the method of manufactured 
solutions. Let 

u(x, y, z) = (y(l - y)z{\ - z), 0, 0), p(x, y, z) = 0.5 - x. 

The right-hand side / is defined accordingly and the corresponding Dirichlet boundary 
conditions are applied via the Nitsche method on the entire boundary Y such that u 



and p solve the Stokes problem ( 1.1 ). 

Let 6 = 0.01 be a perturbation factor. We define three different families of mesh 
configurations, each parametrized over N with h — 1/N, for the background domain 

n*-. 

(A) n* = [-hS, 1 + hS} 3 , divided into N 3 subcubes; 

(B) Q.* = [-h/3, 1 + h/3} 3 , divided into N 3 subcubes; 

(C) fl* = [-h(l - 5), 1 + h(l - S)} 3 , divided into (N + 2) 3 subcubes. 

The final meshes result from tessellating each subcube into 6 tetrahedra. For the 
scenario (A), the background mesh is almost entirely covered by the domain Q; while 
scenario (C) represents the other extreme: the outermost layer of tetrahedra is only 
barely intersected by fi. Scenario (B) illustrates a middle ground. 

For the case V h x Pfc, we take fa = 0.2, fa = 1.0, fa = 0.05 and 7 = 10 as the 



stabilization parameters involved in ( 3.13[ ); while for Vh x P h ' c , we take fa = 0.25 



fa = 0.1 and 7 = 10. To solve the resulting systems of equations, we apply a 
transpose-free quasi-minimal residual (TFQMR) solver with an algebraic multigrid 
preconditioner. The constant pressure mode is filtered out in the iterative solver. We 
observed that the iterative solvers converged in between 6 and 25 iterations. The 
[H 1 (il*)) d error of the velocity approximation and the L 2 (fT) error of the pressure 
approximation were computed, using the natural extensions of the exact solutions to 
f2* , for each mesh configuration and a series o f me sh sizes. 



The resulting errors are plotted in Figure 8.1 and Figure 8.2 for Vh x P,j' dc and 



Vh x Pi, respectively. Theorem 6.1 predicts first order convergence for the iJ 1 -norm 
of the velocity error and the L^-norm of the pressure error. These orders are also 
obtained in the numerical experiments: both for Vh X P^' dc and V/, x P^ and each of 
the three different scenarios, the errors monotonically decrease and seem to converge 
towards zero by (at least) first order. 

We note that the results for the different scenarios illustrate that the positioning 
of background mesh does affect the magnitude of the errors to some extent. For 
the scenario (A), the convergence rates for both the velocity and the pressure seem 
fairly uniform over the range of mesh sizes considered. We observe the same for the 
scenario (B), though the errors and rates are a little higher. For the scenario (C), 
the convergence rates for the L 2 norm of the pressure are somewhat less uniform 
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.1. The case V h X P~ 



errors for the three different mesh configurations (A), (B) and 



(C) versus maximal element diameter h n 



The legend gives the fitted slope for each configuration. 



Top: H -error \\u — n* for the velocity. Bottom: L -error \\p — f or the pressure. 



for the case Vh x P^, and the errors and rates are again higher for both pairs of 
finite clement spaces. As a consequence, we remark that for a series of background 
meshes where the location of the surface varies significantly with respect to the mesh 
configuration, non-monotone decrease of the errors may be observed. We also note 
that superconvergence is observed and is most clearly pronounced in scenario (C). 
This is related to the definition of the norms || • ||i,n* and || • ||q« which extend to 
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the entire fictitious domain f2* . In scenario (C), the fictitious domain fi* extends a 
distance h from the boundary of the computational domain Q. The volume of tt* will 
thus decrease in size during mesh refinement and contribute to the observed rates of 
superconvergence. 
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Scaled condition numbers for Vj, X with varying ghost-penalty stabilization parameters (B = 
h = Ps (each row corresponds to one (5), a varying domain Q = [—Li] 3 and fixed background 
domain Q* = [-1, l] 3 . 
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Table 8.2 

Scaled condition numbers for Vh X P^' dc with varying ghost-penalty stabilization parameter 
j3 = /?2 (each row corresponds to one ft), a varying domain Q = [— £,£] 3 and fixed background 
domain Q* = [— 1, l] 3 . 



8.3. Influence of the boundary position on the condition number. Next, 
we consider a numerical example to demonstrate that the condition number of the 
matrix A corresponding to the stabilized fictitious domain bilinear form, as defined 
by (7.1), is bounded and that the bound is independent of the boundary position 
relative to the background mesh. 

We consider the domain fi* = [— 1, l] 3 tessellated by uniformly dividing the do- 
main into 10 3 cubes, with each cube subdivided into 6 tetrahedra. The domain 
Q = is defined by [—1, 1] 3 , where we have in mind I ranging from 0.9 to 1.0. Note 
that when I is close to 1.0, almost the entire background mesh is included in the com- 
putational domain. On the other hand, as I approaches 0.9, some of the outermost 
elements of the background mesh will only barely intersect Q. So, as I varies between 
1.0 and 0.9, the smallest ratio r of \T P\ Sl| to \T\ for the elements T in the outermost 
layer varies between 1.0 and 0.0. For each I, we compute the condition number of the 
corresponding matrix A, letting /3q = fi\ =0.1, 7 = 10, and varying P2 — P3 — /3. 
The condition number was computed as the ratio of the absolute value of the largest 
(in modulus) eigenvalue and the smallest (in modulus) nonzero eigenvalue of the sym- 
metric matrix A. The resulting condition numbers, scaled by the square mesh size 
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Fig. 8.3. Semilog arithmic plot of the scaled condition number for Vj, X with varying ghost- 
penalty stabilization parameters /3 = fa = fa- 



for a series of j3 and I = 0.99,0.95,0.91,0.901 are given in Table 8.1 and 

,0,dc 



h 2 « 0.35 2 , 

Figure pO|for V h x Pi and in Table^for V h x P°' dc . First, consider the case V h x P\. 
For /3 = 0.0, the scaled condition number is low (386) when I = 0.99; that is, when 
the ratio r is almost 1. However, the scaled condition number increases dramatically 
as I, and hence the ratio r is reduced. Thus, if no ghost-penalty terms are included, 
the scaled condition number seems unbounded as / tends to 0.9. On the other hand, 
in the cases where j3 is positive, the scaled condition number only grows moderately 
as the ratio is significantly reduced and seems bounded. We note however that the 
condition number grows with the penalty param eter for f3 > 0.025. Finally, similar 
observations apply in the case Vh x P°' dc (Table & 



8.4. Stokes flow in a complex geometry. We conclude the section with an 
example of Stokes flow in a computational domain where the boundary is described 
by a complex surface geometry. The geometry is taken from a part of an arterial 
network known as the Circle of Willis which is located close to the human brain. It is 
known that the network is prone to develop aneurysms and therefore the computer- 
assisted study of the blood flow in the Circle of Willis has been a recent subject 
of interest, see for instance Isaksen et al. [21], Steinman et al. |34j . Valen-Sendstad 
et al. [35] - However, the purpose of this example is not to perform a realistic study 
of the blood flow dynamics. Rather, we would like to demonstrate the principal 
applicability of the developed method to simulation scenarios where complex three- 
dimensional geometries are involved. The extension of the work to numerically solve 
the time-dependent Navier-Stokes equations in a biomedical relevant regime is the 
subject of future research. 

The blood vessel geometry is embedded in a structured background mesh as 
illustrated in Figure 8.4 As before, the velocity is prescribed on the entire boundary 
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r where we set it = on the arterial walls and u = 1200 mm/s on the inlet boundary. 
The two outflow velocities were set in such a way that total flux was balanced. 

The pressure and velocity approximation as computed on the fictitious domain 



mesh T* are shown in Figure 8.4 and 8.5 respectively. Although the fictitious domain 
mesh T* provides only a coarse resolution of the aneurysm geometry, the values of 
the velocity approximation clearly conforms to the required boundary values on the 
actual surface geometry. 

9. Conclusions. We have presented a stabilized finite element method for the 
solution of the Stokes problem on fictitious domains and proved optimal order con- 
vergence. The theoretical convergence rates have been verified numerically. We have 
also proved that the condition number of the stiffness matrix remains bounded, inde- 
pendently of the position of the fictitious boundary relative to the background mesh. 

While we have here restricted our attention to the static Stokes model problem, 
the main motivation for the methodology and implementation presented in this paper 
is for the treatment of the time-dependent Navier-Stokes equations and, ultimately, 
fluid-structure interaction on complex and evolving geometries. We address this issue 
in future work. 
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